## Heterogeneity Analysis ##

pool_state <- pool %>% mutate(
  state = case_when(
    state == 9 ~ "Uttar Pradesh",
    state == 10 ~ "Bihar",
    state == 19 ~ "West Bengal",
    state == 20 ~ "Jharkhand",
    state == 21 ~ "Odisha",
    state == 23 ~ "Madhya Pradesh"
  )) %>% filter(state == "Bihar")

# With Village Fixed Effects ---------------------------------------------------

# LPG Connection (Binary)
model.dd.vfe.uselpg.scst = fixest::feols(uselpg ~ scst*wave + age + gender +
                                           religion + education + hhsize +
                                           lnexp + bplaay | village, 
                                         data = pool_state, weights = pool_state$weights)

# LPG Home Delivery (Binary)
model.dd.vfe.convenience.scst = fixest::feols(convenience ~ scst*wave + age + gender +
                                                religion + education + hhsize +
                                                lnexp + bplaay | village, 
                                              data = pool_state, weights = pool_state$weights)

# Grid Connection (Binary)
model.dd.vfe.usegrid.scst = fixest::feols(usegrid ~ scst*wave + age + gender +
                                            religion + education + hhsize +
                                            lnexp + bplaay | village, 
                                          data = pool_state, weights = pool_state$weights)

# With Household Fixed Effects -------------------------------------------------

# Regression Models for Energy Access Response Variables
# Robust standard errors are clustered at household level.

# LPG Connection (Binary)
model.dd.hfe.uselpg.scst = fixest::feols(uselpg ~ scst*wave + hhsize +
                                           lnexp + bplaay | hh, 
                                         data = pool_state, weights = pool_state$weights)

# LPG Home Delivery (Binary)
model.dd.hfe.convenience.scst = fixest::feols(convenience ~ scst*wave + hhsize +
                                                lnexp + bplaay | hh, 
                                              data = pool_state, weights = pool_state$weights)

# Grid Connection (Binary)
model.dd.hfe.usegrid.scst = fixest::feols(usegrid ~ scst*wave +  hhsize +
                                            lnexp + bplaay | hh, 
                                          data = pool_state, weights = pool_state$weights)

## Model and coefficient lists for household fixed effects combined ##
models.hfe.dd_access_scst <- list(model.dd.hfe.uselpg.scst, 
                                  model.dd.hfe.convenience.scst, 
                                  model.dd.hfe.usegrid.scst)

## Model and coefficient lists for village and household fixed effects combined ##
models.vfe.hfe.dd_access_scst <- list(model.dd.vfe.uselpg.scst, model.dd.vfe.convenience.scst,
                                      model.dd.vfe.usegrid.scst, model.dd.hfe.uselpg.scst, 
                                      model.dd.hfe.convenience.scst, model.dd.hfe.usegrid.scst)

esttex(models.vfe.hfe.dd_access_scst, 
       se = "cluster", 
       digits = 3,
       dict = c(scst = "SC/ST", wave = "2018", uselpg = "LPG Connection",
                usegrid = "Grid Connection", convenience = "LPG Home Delivery",
                village = "Village", hh = "Household"),
       keep = c(":", "2018"),
       fitstat = c("r2"),
       fixef_sizes = T,
       float = F,
       replace = T,
       file = c("./Manuscript/Tables/dd_access_scst_state_BI.tex"))

# State-wise analysis results collection ---------------------------------------

state_bi_fits <- models.hfe.dd_access_scst
